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Abstract 

Betula platyphylla Suk (birch) is a fast-growing woody species that is important in pulp industries and the biofuels. However, 
as an important pulp species, few studies had been performed on its wood formation. In the present study, we investigated 
the molecular responses of birch xylem to artificial bending and gravitational stimuli. After trunks of birch trees were 
subjected to bending for 8 weeks, the cellulose content was significantly greater in tension wood (TW) than in opposite 
wood (OW) or normal wood (NW), whereas the lignin content in TW was significantly lower than that in OW and NW. In 
addition, TW grew more rapidly than OW and generated TW-specific fibers with an additional G-layer. Three transcriptome 
libraries were constructed from TW, OW and NW of S. platyphylla, respectively, after the plants were subjected to artificial 
bending. Overall, 80,909 nonredundant unigenes with a mean size of 768 nt were assembled. Expression profiles were 
generated, and 9,684 genes were found to be significantly differentially expressed among the TW, OW and NW libraries. 
These included genes involved in secondary cell wall structure, wood composition, and cellulose or lignin biosynthesis. Our 
study showed that during TW formation, genes involved in cellulose synthesis were induced, while the expression of lignin 
synthesis-related genes decreased, resulting in increased cellulose content and decreased lignin levels in TW. In addition, 
fasciclin-like arabinogalactan proteins play important role in TW formation. These findings may provide important insights 
into wood formation at the molecular level. 
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Introduction 

Trees are subjected to various constraints during development, 
including environmental constraints such as soil instability and 
damage caused by the actions of wind or snow, as well as 
biological constraints (such as limited light), which induce 
gravitropic and phototropic responses. To remain upright under 
adverse conditions, angiosperms have developed specialized wood 
tissue called tension wood (TW), which is formed in response to 
gravity. TW is formed on the upper sides of branches and leaning 
trunks of angiosperms and exerts a contractile force that allows a 
tree trunk to remain vertical or a branch to remain horizontal [1]. 
The other type of reaction wood is called compression wood (CW), 
which is formed on the ground-facing sides of trunks and the 
leaning stems of gymnosperns [2]. There are two major methods 
for artificially inducing the formation of tension wood, including 
bending the stem into a loop (extreme induction) and bending a 
stem using an immobile support. Both TW and opposite wood 
(OW) are generated with gravitropic stimuli, mechanical stimuli or 
both [2] . Compared with OW or normal wood (NW), TW exhibits 
some unusual characteristics [3]. Secondary cell walls of TW 
contain a gelatinous layer (G-layer), which consists of crystalline 
cellulose at a low microfibrillar angle [4] . The presence of a G- 
layer in TW results in the higher cellulose content and lower 



proportion of lignin found in TW compared with OW or NW. In 
addition, the G-layer also has a reduced proportion of vessel 
elements, and it is often formed more quickly than NW [5]. The 
formation of OW is either strongly retarded or completely 
inhibited. 

To date, few studies have addressed the chemical composition 
of OW [6]. Some studies have demonstrated that OW has the 
same chemical composition as NW [7]. A comparison of gene 
expression during the formation of TW and OW would help to 
elucidate the process of wood formation. Andersson-Gunneras 
et al. generated an expressed sequence tag library from TW- 
forming tissues in Populus tremula (L.) tremuloides (Michx.) and 
studied gene expression profiles during TW formation in Populus 
tremula (L.) in two growing seasons using microarray and 
metabolite analyses [8] . Jin et al. constructed a cDNA library 
from secondary xylem in the stem of a 2-year-old yellow poplar 
that was subjected to bending for 6 h at a 45° angle [9]. In total, 
5,982 high-quality expressed sequence tags (ESTs) were generated, 
which were clustered into 1,733 unigenes, including 822 contigs 
and 911 singletons. Analysis of these ESTs identified many genes 
that are involved in cell wall biosynthesis and modification. Qju 
et al. studied gene expression in Eucalyptus nitens branches bent at a 
45° angle using microarrays that contained 4,900 cDNAs from 
xylem, and important genes involved in responses to gravitational 
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stxess in eucalyptus xylem were identified [10]. To characterize 
gene expression during TW formation in eucalyptus, cDNA array 
technology was employed. In total, 196 genes were found to be 
differentially regulated between control and bent wood samples. 
Some of these genes displayed distinctive expression patterns 
related to changes in secondary cell wall structure and composi- 
tion. Analysis of these gene expression profiles provided new 
insights into the regulatory network of genes that are preferentially 
expressed in xylem [II]. Although these studies identified genes 
that are differentially expressed in TW versus controls, only a 
limited number of genes were investigated. However, little is 
known about the expression of genes associated with cell wall 
biosynthesis and modification on a large scale. 

Next-generation (NG) sequencing techniques are low cost, high 
throughput sequencing methods that can generate information 
about numerous expressed genes in a short period of time. 
Platforms for NG sequencing include the Genome Analyzer 
(Solexa/IUumina), 454 (Roche) and ABI-SOLiD (Applied Biosys- 
tems) [12]. RNA-Seq is an effective approach to transcriptome 
profiling that uses NG sequencing technologies. RNA-Seq 
provides a far more precise measurement of the transcript levels 
of genes and their isoforms than other methods of analysis on a 
genomic scale [13]. Morco\-cr, RN.\-Seq can provide absolute 
gene expression measurements rather than relative measurements, 
thus overcoming many limitations of microarray analysis [14,15]. 

Betula platyphylla Suk (birch) is a fast-growing woody species that 
it is tolerant to Ught, drought and flooding and adapts well to 
numerous t)pes of soil. This widely grown tree is important in the 
biofuels and pulp industries. However, few essential regulatory 
genes have been identified to date. Such genes could be used for 
the genetic improvement of birch in particular to enhance its value 
as a source of biofuels and biomaterials. The aim of the present 
study was to investigate the molecular responses of birch to tension 
stimulus and lay the foundations for improving the wood 
properties of birch planted for the production of biofuels and 
pulp. The ceUulose and lignin contents and wood anatomy of TW, 
OW and NW were analyzed in birch stems subjected to 8 weeks of 
mechanical bending. The results suggest that tension stimulus can 
influence xylem development in birch. As genes that are expressed 
early in the formation of reaction wood may significantly influence 
the subsequent quality of the wood, we examined the transcrip- 
tomes of TW, OW and NW in B. plalyphylla subjected to bending 
for 2 weeks. We determined a complete set of birch xylem genes 
that were differentially expressed in response to gravitational stress 
and mechanical bending, and genes were identified that are 
involved in the biosynthesis of wood components, such as cellulose 
and lignin. The results of our study are potentially useful for 
improving wood properties in plants via genetic engineering. 

Materials and Methods 

Ethics Statement 

The location for growing the birch trees used in this study is not 
required specific permission, because it is the experimental field of 
our Laboratory (State Key Laboratory of Tree Genetics and 
Breeding, Northeast Forestry University), which is named as 
"orchard of birch intensive breeding", located in Northeast 
Forestry University. In this orchard, we have the experimental 
field, and all the birch trees used in this study were planted by us; 
therefore, we have the right to use these birch trees in study. 

Plant Material 

Beginning in May 2011, 2-year-old birch trees were grown 
under natural conditions in a field in Harbin (the People's 



Republic of China), located at 45°44'N and 126°36'E. All of the 
trees were well watered. TW formation was induced in trees 
during the most active period of growth by bending the tree and 
tying it to its neighbor tree in such a way that the midpoint of the 
stem was at an angle of approximately 45°; upright trees were used 
as the control. Two weeks after applying the stimulus, xylem 
samples were collected from two bent trees and two control trees. 
Tissues were collected by peeling the bark and scraping the 
exposed difiFerentiating xylem (2—3 mm in depth) from both the 
upper (TW) and lower (OW) sides of the bent stems with a scalpel, 
as well as from corresponding locations on the controls (nonbent 
trees, NW). Xylem samples were immediately frozen in liquid 
nitrogen and stored at — 80°C. 

Chemical Content and Wood Anatomy Analysis 

Wood samples treated as described above were also harvested 
after 8 wks of bending to produce material for transverse sections 
(10-15 |J,m thick) and to analyze the cellulose and hgnin contents. 
To characterize the presence of TW by the appearance of G- 
layers, sections were subjected to environmental scanning electron 
microscopy (ESEM). The cellulose and lignin contents were 
determined with a Fibertec^^' 2010 & M6 system (Foss, Sweden) 
following the manufacturer's protocol (www.foss.dk). Analyses 
were performed in triplicate. Analysis of variance (ANOVA) was 
performed using SPSS software. 

Construction of lllumina Library and Solexa Sequencing 

Total RNA was isolated from OW, TW and NW using the 
CTAB method and digested with DNase I (RNase free) to remove 
DNA contamination. Since each treatment sample contains two 
trees, RNA was respectively isolated from these two trees and then 
pooled for library construction. lllumina Solexa sequencing using 
the GAII platform was carried out at the Beijing Genomics 
Institute (BGI; Shenzhen, China). Briefly, mRNA was isolated 
from total RNA from each sampk- using poly-T oligo-attached 
magnetic beads (lllumina) and sheared into small fragments using 
divalent cations (lllumina) at 94°C for 5 min. Complementary 
DNA (cDNA) was synthesized using mRNA fragments as 
templates and random primers as reverse complementary primers. 
To create blunt ends, the synthesized cDNAs were digested with 
T4 DNA polymerase and Klenow DNA polymerase. The 3' ends 
of the blunt phosphorylated cDNA fragments were added with an 
adenosine (A) base using the polymerase activity of the Klenow 
fragment. Then, specific adaptors equipped with single 'T' base 
overhangs at their 3' ends were ligated to the cDNA fragments 
supplemented with the 'A' base. The ligation products were run on 
a gel to select fragments of approximately 200 bp in length. These 
selected fragments were PGR amplified using a primer set that 
specifically annealed to the ends of the adapters. The PGR 
products were then purified with a QIAquick PGR Purification 
Kit that was used for construction of the lllumina library. After 
evaluation of the size, purity and contx^ntration of the constructed 
library, single-stranded cDNA fragments were annealed to the 
flow cell surface in a cluster station (lllumina, San Diego, CA). 
Sequencing-by-synthesis was performed with an lllumina Genome 
Analyser, and the reads that passed default quality filtering 
parameters in the lllumina GA Pipeline GERALD stage were 
retained for further analysis. 

Data Analysis 

High quahty reads were generated by removing adaptor 
sequences, reads with only adaptor sequences, reads with 
unknown sequences 'N' and reads with a copy number of one 
(probably sequencing error). Transcriptome de novo assembly was 
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performed in each library using the short oligonucleotide analysis 
program SOAPdenovo [16]. SOAPdenovo combines short reads 
with a certain length of overlap to form contigs, and the reads are 
mapped back to the contigs. From the paired-end reads, contigs 
from the same transcript and the distances between these contigs 
were both detected. Scaffolds were formed by connecting the 
contigs from the same transcript using SOAPdenovo, with N 
representing unknown sequences between the two contigs. Paired- 
end reads were employed for gap fdling within scaffolds to 
generate sequences with the fewest Ns that could not be extended 
at either end. Such sequences were considered to be tentative 
unigenes (TUGs). The TUGs from each library were further 
assembled from a single set of nonredundant unigenes using 
TGICL software [17]. The abundance of each gene was 
determined by calculating the reads per kb per million reads 
(RPKM) value, as described by Mortazavi et at. [18]. 

Functional Annotation 

Nonredundant unigenes were subjected to BLASTX analysis 
against the NR of NCBI and Swiss-Prot protein databases to 
search for similarity. The unigenes with BLASTX E-value > 1 0 ^ 
were discarded during functional annotation. Gene Ontology 
(GO) annotation was assigned using the Blast2go program [19]. 
After GO annotation was assigned for each unigene, WEGO 
software was used to perform GO functional classification for all 
unigenes according to molecular function, biological process and 
cellular component [20]. 

Real-time Quantitative RT-PCR Analysis 

To test the overall validity of the expression levels of genes 
determined by DGE analysis, real time RT-PCR was performed. 
One microgram of total RNA from each sample was reverse 
transcribed using ohgo-deoxythymidine as a primer in a 10 [ll 
volume. Real time RT-PCR was performed using an MJ 
Opticon''"'^ ^ machine (Bio-Rad, Hercules, CA). Actin (GenBank 
number: HO 112155) and ubiquitin (GenBank number: HO 112156) 
genes were used as internal controls. The reaction mixture (25 |J.l) 
contained 12.5 |xl of SYBR Green Real-time PGR Master Mix 
(Toyobo), 0.5 |J.M of forward or reverse primers (Table 1 in File 
SI), and cDNA template (equivalent to ICQ ng of total RNA). The 
amplification procedure was performed using the following 
parameters: 94°C for 30 s; 45 cycles at 94°C for 12 s, 58°C for 
30 s, 72°C for 45 s and 80°C for 1 s for plate reading. A melting 
curve was generated for each sample to assess the purity of the 
amplified products. AU experiments were carried out with three 
biological replicates. The expression levels were calculated from 
the threshold cycle using the delta-delta CT method. 

Results and Discussion 

Cellulose and Lignin Content and Wood Anatomy 

The cellulose and lignin contents of TW, OW and NW were 
analyzed in the stems of 2-year-old birch trees subjected to 
mechanical bending for 2 and 8 weeks. There was no significant 
difference in the cellulose and lignin contents among the three 
tissues subjected to bending for 2 weeks (data not shown). 
However, after bending for 8 weeks, the cellulose content was 
significantly higher in TW than in OW or NW, whereas the lignin 
content was significantly lower in TW than in the other types of 
wood (Figure 1). 

Samples were collected from the upper (TW) and lower (OW) 
sides of bent 2-year-old trees and from straight trees (NW) after 8 
weeks of bending. Sections were produced from these samples and 
examined. There were obvious differences in growth rates between 
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Figure 1. Cellulose and lignin contents of TW, OW and NW. * 

Significant differences (P<0.05) between TW and OW or NW. Cellulose 
and lignin content are expressed as a percentage of cell wall residue 
(CWR). 

doi:10.1371/journal.pone.0087566.g001 

TW and OW, which led to the unusual growth of birch stems 
(Figwe 2M, N). In addition, TW exhibited characteristic fibers 
with an additional G-layer (Figure 2A-D) and contained fewer 
vessels with smaller lumens than OW (Figure 2E-H) or NW 
(Figwe 2I-L). No G-layer was present in NW, or in OW subjected 
to bending for 8 weeks (Figure 2E-L). The increased cellulose 
content and fiber numbers in response to tension is thought to be 
part of a stem-sensing mechanism, which leads to increased tensile 
strength to counteract the gravitational forces acting on the 
leaning stem. 

De novo Transcriptome Assembly 

The gene expression in response to artificial bending were 
usually studied after bending for 6 h to 3 weeks [8,9,1 1], and that 
period is important for study the tension wood formation. Three 
cDNA libraries were constructed from B. platyphylla TW, OW and 
NW, after bending for two weeks (referred to as the TW, OW and 
NW libraries). In total, 52,139,846, 54,878,244 and 54,696,478 
clean reads were generated from the TW, OW and NW libraries, 
respectively. The data of clean reads can be accessed at the NCBI 
Sequence Read Archive under the accession number SRA053683 
(SRR5 13429, SRR520271, SRR520272). The reads from each 
library were assembled into contigs using Trinity software [21], 
and the contigs were further assembled into tentative unigenes 
(TUGs). In total, 77,783, 66,909 and 74,068 TUGs were 
generated from the TW, OW and NW libraries, respectively 
(Table 2 in File SI). Finally, all unigenes from the three libraries 
were assembled together to form a single set of non-redundant 
unigenes using TGICL software [17]. A total of 80,909 non- 
redundant unigenes (GenBank No.: KA198535-KA279443) was 
generated, with a mean size of 768 nt and a total length of 
62,112,335 nt. The distribution of unigene lengths is shown in 
Table 1. 

Annotation of Nonredundant Unigenes 

Among the non-redundant genes, 43,245 and 29,977 genes 
showed significant similarity to known proteins in the Nr and 
Swiss-Prot databases, respectively. 

GO classification was performed to study the functional group 
distribution of non-redundant unigenes (Table 3 in File SI & 
Figure SI). In all three libraries, cellular components were the 
most abundant category, followed by biological processes and 
molecular function. Of the three principal gene ontologies, the 
most abundant subcategories in the three libraries were 'all cell', 
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Figure 2. Images of wood anatomy of TW, OW and NW. Environmental scanning electron microphotographs of the transverse surface of 
tension wood (A-D) and opposite wood (E-H) in stems (LS) of birch cultivated for 8 wks under control conditions (l-L). A, E, I: x300; B, F, J: x500; C, G, 
K: xlOOO; D, H, L: x3000. Scale bars are described in O. V, vessel; F, fibre; GF, G-layer fibre. IVl and N show the unusual growth of different parts of 
birch stems subjected to bending for 8 wks. 
doi:1 0.1 371 /journal.pone.0087566.g002 



Table 1. Length distribution of all unigenes. 



Length of non-redundant 
unigenes 


Total Number 


Percentage 


100-500 nt 


45,950 


56.79% 


500-1000 nt 


15,524 


19.19% 


1000-1500 nt 


7,945 


9.82% 


1 500-2000 nt 


4,940 


6.11% 


> = 2000 nt 


6,550 


8.09% 


N50 [1] 


1309 




Mean 


768 




All Unigene 


80,909 




Length of all Unigene (nt) 


62,112,335 





*N50 = median length of all unigenes. 
*Mean = average length of all unigenes. 
doi:l 0.1 371/journal.pone.0087566.t001 



'cell part', 'catalytic activity', 'binding', 'organelle', 'metabolic 
process' and 'cellular process', suggesting that these functional 
processes play major roles in TW, OW and NW. 

The Most Abundant Genes 

Genes that are abundant in certain tissues or development 
stages may be important for the tissue or developmental stage; 
therefore, we identified the ten most abundant genes in the TW, 
OW and NW libraries (Table 2). Seven transcripts were highly 
abundant in all of the libraries, including two lipid transfer 
proteins (LTP), allergenic isoflavone reductase-like, agl3, poly- 
ubiquitin, cytochrome P450 like_TBP and hypothetical protein 
(Table 2). These findings suggest that these genes play essential 
roles in wood formation. Genes encoding LTP, cytochrome P450 
like_TBP, rRNA intron-encoded homing endonuclease, phi-1 and 
C3HL domain class transcription factor were abundant in OW, 
with much less abundance in TW, suggesting that these genes play 
important roles in OW development rather than TW develop- 
ment. On the contrary, genes encoding FLA, glycine-proline rich 
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Table 2. Expression profiles of the top 10 abundant genes in TW, OW and NW. 





GenBank Number: 


RPKM 






Functional annotation 




TW 


OW 


NW 




NW 


KA275752 


3724.332 


15053.39 


1 3608.72 


lipid transfer protein family protein 


KA250254 


7559.774 


9818.979 


11268.71 


allergenic isoflavone reductase-like 


KA244244 


6013.817 


11912.6 


9884.284 


lipid transfer protein (LTP) family protein 


KA200383 


8000.385 


6803.264 


6643.729 


agl3 


KA262474 


2025.879 


2844.252 


5724.919 


xylem sap protein 10 kDa 


KA201101 


5649.985 


4611.468 


5311.726 


polyubiquitin 


KA259905 


3619.965 


13164.14 


5246.218 


cytochrome P450 like_TBP 


KA250879 


4966.727 


5211.61 


4390.03 


hypothetical protein 


KA257282 


6723.616 


2745.166 


3908.591 


glycine-proline rich protein 


KA257879 


2550.106 


2733.757 


3309.829 


polyubiquitin containing 7 ubiquitin monomers 


TW 


KA200383 


8000.385 


6803.264 


6643.729 


agl3 


KA250254 


7559.774 


9818.979 


11268.71 


allergenic isoflavone reductase-like protein 


KA219691 


7317.225 


242.5356 


654.6115 


fasciclin-like arabinogalactan protein 


KA257282 


6723.616 


2745.166 


3908.591 


glycine-proline rich protein 


KA244244 


6013.817 


11912.6 


9884.284 


lipid transfer protein (LTP) family protein 


KA201101 


5649.985 


4611.468 


5311.726 


polyubiquitin 


KA250631 


5167.007 


2003.051 


2410.727 


class II chitinase 


KA250879 


4966.727 


5211.61 


4390.03 


hypothetical protein 


KA275752 


3724.332 


15053.39 


13608.72 


lipid transfer protein family protein 


KA259905 


3619.965 


13164.14 


5246.218 


cytochrome P450 like_TBP 


OW 


KA275752 


3724.332 


15053.39 


13608.72 


lipid transfer protein family protein 


KA259905 


3619.965 


13164.14 


5246.218 


cytochrome P450 like_TBP 


KA244244 


6013.817 


11912.6 


9884.284 


lipid transfer protein (LTP) family protein 


KA250254 


7559.774 


9818.979 


11268.71 


allergenic isoflavone reductase-like protein 


KA259904 


2660.703 


7896.881 


3299.789 


rRNA intron-encoded homing endonuclease 


KA200383 


8000.385 


6803.264 


6643.729 


agl3 


KA250879 


4966.727 


5211.61 


4390.03 


hypothetical protein 


KA201101 


5649.985 


4611.468 


5311.726 


polyubiquitin 


KA249327 


954.0103 


2977.693 


845.0917 


phi-1 


KA262532 


779.1426 


1631.054 


1012.204 


C3HL domain class transcription factor 



doi:1 0.1 371 /journal.pone.0087566.t002 



protein and class II chitinase were higlily abundant in TW, with 
little expression in OW; obviously, these genes play important 
roles in TW formation. 

Identifying Differentially Expressed Genes in TW, OW and 
NW 

Unigenes more than 400 bp in length, FDR <0.0001 and the 
absolute value of |log2Ratio| &2, were identified as significandy 
differentially expressed using DGE analysis (Table SI). The 
distributions of differentially expressed genes among the three 
libraries are shown in Table 3. GO classification was performed to 
study the functional group distribution of the differentially 
expressed genes (Figure 3). Among the three principal gene 
ontologies, the most abundant subcategories in TW, NW and OW 
were ^aU cell', 'cell part', 'catalytic activity', 'binding', 'organelle', 
'metabolic process' and 'cellular process', thereby suggesting that 



these functional processes play major roles in the formation of TW 
and OW. 

Fasciclin-like Arabinogalactan Proteins are Abundant and 
Upregulated in TW 

Our results show that a fasciclin-like arabinogalactan protein 
(FLA) gene (KA2 19691) is one of the top 10 most abundant genes 
in the TW library and is strongly induced in TW relative to OW 
and NW (Table 2). In addition, 19 unique FLAs genes were 
identified, and 10 FLAs were induced in TW compared with NW 
and OW (Table 4 in File SI & Figure 4A), suggesting that these 
genes are involved in TW formation. FLAs are a subclass of AGPs 
containing domain(s) that are involved in wood formation [22] and 
the building of the G-layer [23]. Taken together, these results 
suggest that these FLAs are involved in building the G-layer in TW 
of birch. 
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■ T vs. O no vs. N ax vs. N 
molecular function 



antioxidant activity 
binding 
catalytic activity 
emyme regulator activity Q 
molecular transducer activity 
protein binding transcription factor activity 
translation regulator activity 
transporter activity 



O 500 lOOO 1500 2000 2500 3000 



cellular component 



cell 

cell junction 
cell part 
extracellular region 
extracellular region part 
macromolecular complex 
membrane-enclosed lumen 
organelle 



organelle part [ , I ^ 



O 5()0 I (KK) 1 50O 2000 2500 3000 3500 



biological process 

biological adhesion 



biological regulation [.. '■ 
cellular component organization or biogenes is \ \ ' 
cellular process [ '■ 
death Q 
developmental process [ | ' 
establishment of localization \ _ 
growth Q 
immune system process | 

localization \ '■ 
locomotion 
metabolic process [ '■ 
multicellular organismal process |^^^ ' 
multi-organism process ^ 
negative regulation of biological process ^ 
pigmentation 
positive regulation of biological process ^ 
regulation of biological process 
re production 
reproductive process 
response to stimulus |— 
rtiythmic process 
signaUng 
viral reproduction 



O 



Figure 3. Gene ontology (GO) categories assigned to differentially expressed genes in response to artificial bending. The x-axis 
shows the number of unigenes, and the y-axis shows the GO subcategories. 
doi:1 0.1 371 /journal.pone.0087566.g003 
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Table 3. The distributions of differentially expressed genes. 



Number of genes 







TW vs. 


OW vs. 


OW vs. 






NW 


NW 


TW 


Total 




6264 


4533 


4025 


Up-regulated 




2792 


2322 


2470 


Expressed only In 


this tissue 


525 


313 


323 


Down-regulated 




3472 


2211 


1555 


Expressed only in 


control tissue 


501 


315 


480 



*The control tissue refers to the latter tissue in every column. 
doi:l 0.1 371/journal.pone.0087566.t003 

LTP may Play Specific Roles in TW Formation 

LTPs (lipid transfer proteins) are involved in promoting plant 
cell wall loosening and thus may play a role in cell expansion and 
plant growth [24]. We identified 10 unique LTP genes, and 2 LTP 
genes (KA275752 and KA244244) were among the top 10 most 
abundant genes in all three libraries (Table 2), suggesting that they 
play crucial roles in the formation of TW, OW and NW. In 
addition, 7 LTP genes had lower levels of expression in TW than 
in OW or NW (Table 4 in File S1& Figure 4B). Combined with 
the smaller lumen in TW cells, the decreased expression of these 
LTPs in TW may contribute to decreased cell wall loosening or 
inhibited cell expansion during TW formation. 

Genes Involved in Cellulose Synthesis Were Induced 
during TW Formation 

Cellulose, which comprises long chains of (5-1, 4-linked glucosyl 
residues, is formed from uridine diphospho (UDP)-glucose 
catalyzed by cellulose synthase (CesA, EC2.1.4.12) [25]. Our 
results show that 12 CesA genes were highly expressed in TW 
relative to in OW or NW (Table 4 in File SI, Figure 4C & Figure 
S2a). The unigenes KA240937, KA240936 and KA207874 shared 
high identities with the irregular xylem3 (irx3) from Arabidopsis, 
which was required for secondary cell wall synthesis [26], and all 
of the three genes were upregulated in TW compared those in 
OW or NW (Table 4 in File SI & Figure 4C). Moreover, our 
results show that cellulose levels in TW are elevated compared 
with in OW and NW (Figure 1). Taken together, these results 
suggest that the increased expression level of CesA results in 
elevated cellulose levels; therefore, these CesA genes play roles in 
the synthesis of cellulose and are involved in TW formation. 

Cellulose biosynthesis rosettes use UDP-glucose as a substrate 
[27], and the UDP-glucose synthesis pathways involve a serial 
enzymes, such as sucrose synthase (SuSy), sucrose phosphate 
synthase (SPS), UDP-glucose 4-epimerase (GALE), glucose- 1- 
phosphate uridylyltransferase and a-expansins [28-32]. Our 
results show that these genes were upregulated in TW than in 
NW or OW (Table 4 in File SI, Figure 4C & Figure S2a). For 
instance, one SuSy genes was significandy induced, and three SuSji 
were slightly enriched in TW compare with OW or NW (Table 4 
in File SI & Figure 4C). Six SPSs were identified, and 5 of them 
were highly induced in TW compared with OW or NW. Two 
GALE genes were highly enriched in TW. Two glucose- Tphosphate 
uridylyltransferase genes were significandy enriched in TW compared 
with in OW or NW. Moreover, 3 a-expansin genes (KA276137, 
KA202593 and KA2 70930) were highly induced in TW relative to 
in OW (Table 4 in File SI & Figure 4C). The fact that many 
cellulose biosynthesis-related genes were induced in TW, com- 



bined with the high levels of cellulose that accumulated in TW, 
clearly suggest that cellulose biosynthesis is highly activated in 
TW. Therefore, the maintenance of high levels of cellulose 
biosynthesis is important for TW formation. 

Genes involved in the Shikimate, Phenylpropanoid and 
Monolignol Biosynthesis Pathways were Downregulated 
in TW 

The lignin content decreased significantly in TW subjected to 
bending for 8 weeks compared with NW and OW (Figure 1), 
which indicates that the lignin biosynthesis pathway is inhibited in 
TW during the bending stress period. The shikimate pathway 
functions in lignin biosynthesis by channeling the flow of carbon 
from sugar metabolism to the biosynthesis of phenylalanine. 
Phenylalanine, the precursor of monolignol, is converted into 
monolignol via the phenylpropanoid and monolignol biosynthesis 
pathways [33,34]. The shikimate pathway includes the enzymes 
dehydroquinate dehydratase-shikimate dehydrogenase (DHQ; 
SDH), shikimate kinase (SK), 5-enolpyruvylshikimate-3-phosphate 
synthase (EPSPS) and others. Most of genes involved in the 
shikimate pathway were downregulated in TW-forming tissues 
relative to OW and NW (Table 4 in File SI, Figure 4D & Figure 
S2b). For instance, 7 DHQ^SDH genes were downregulated, while 
only 2 DHQjSDH genes were upregulated in TW relative to OW 
or NW. Two SK genes were significandy downregulated in TW 
compared with in OW or NW. Two EPSPS genes were 
downregulated in TW compared with OW or NW. In the TW 
of birch, the expression of DHQjSDH, SK and EPSP were all 
downregulated, and lignin content in TW was also decreased 
(Table 4 in File SI, Figure 4D & Figure 1). Taken together, these 
results suggested that the shikimate pathway is positively regulated 
during lignin biosynthesis in both angiosperm and gymnosperm 
plants. 

The biosynthesis of monolignols proceeds through the phenyl- 
propanoid pathway, which starting with the domination of 
phenylalanine to produce cinnamic acid, which then forms 
monolignols via the monolignol biosynthetic pathway. The genes 
involved in the early phenylpropanoid pathway were downregu- 
lated in TW compared with in OW and NW (Table 4 in File SI, 
Figure 4D & Figure S2b). For instance, PAL (phenylalanine 
ammonialyase) is the first enzyme in the phenylpropanoid 
pathway, and four PAL genes (KA266994, KA273229, 
KA277908 and KA243885) were expressed with lower expression 
levels in TW than in OW or NW. The enzyme 4CL (4-coumarate- 
coa ligase ) is the key enzyme in the phenylpropanoid pathway, 
and eight 4CLs were downregulated in TW compared with in NW 
or OW. In addition three CCoAOMT (caflFeoyl-CoA O-methyl- 
transferase) and six CAD (cinnamyl alcohol dehydrogenase) genes 
were highly downregulated in TW relative to in NW or OW. The 
decreased transcript levels of genes involved in the phenylpropa- 
noid and monolignols biosynthesis pathways combined with 
decrease of lignin content suggest that lignin biosynthesis is 
inhibited during TW formation, and these genes may play 
important roles during TW formation in birch. 

POD (peroxidase) and laccase are involved in monolignol 
polymerization and modifications. The transcript levels of 20 POD 
genes were downregulated in TW compared with in NW or OW. 
The transcript levels of 4 laccase genes (KA265956, KA201153, 
KA207714 and KA211912) were significandy decreased in TW 
compared with in OW or NW. These results indicate that 
monolignol polymerization and modifications are inhibited during 
TW formation. 

From the above results, it is evident that the genes involved in 
the shikimate, phenylpropanoid and monolignol biosynthesis 
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Figure 4. Hierarchical cluster analysis of different gene families with altered expression levels in response to artificial bending. All 

ratios are log2 transformed so that inductions and repressions of identical magnitudes are numerically equal but opposite in direction. Log ratios of 0 
(ratios of 1) are shown in black, and increasingly positive (induction) or negative (repression) log ratios are shown in red or green, respectively, with 
increasing intensities. Red indicates induction, and green indicates repression in arrays. A: FLA family genes, B: LTP family genes, C: Cellulose 
biosynthesis related genes, D: Lignin biosynthesis related genes, E: TFs, F: Growth hormone-related genes. 
doi:1 0.1 371 /journal.pone.0087566.g004 



pathways are generally downregulated in TW. In addition, the 
lignin content was decreased in the TW of birch. Taken together, 
these results suggest that the lignin biosynthesis is inhibited in TW 
by the downregulation of genes involved in the shikimate, 
phenylpropanoid and monolignol biosynthesis pathways. 

Identification of Wood-related Transcription Factors 

We focused our study on NAC and MYB TFs, because they 
play important roles in secondary cell wall modification [35-38]. 
Using hierarchical cluster analysis, differentially expressed JVAC 
genes were assigned into two groups according to their expression 



patterns (Table 4 in File SI & Figure 4E): one group contained the 
MACs that were upregulated in TW and another included NACs 
that were downregulated in TW compared with OW or NW. The 
Unigenes KA246202 and KA269898 matched VMDl from 
Arahidopsis, and both of them were upregulated in TW compared 
in OW or NW (Table 4 in File SI & Figure 4E). The Unigene 
KA255359 matched Arahidopsis NSTl, a key regulator of the 
formation of secondary walls [39], and its transcript level was 
altered significantly between OW, TW and NW, suggesting that it 
is involved in the formation of secondary walls. The changes in 
NAC expression among TW, OW and NW suggest that NACs are 
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closely associated with TW formation, and they may play roles in 
mediating TW development. MTB family proteins play positive or 
opposite roles in the lignin biosynthetic process [40,41]. Nearly- 
half of the A4TB.S were upregulated in TW vs. OW or NW, and the 
remaining MYBs were downregulated in TW compared with OW 
and NW (Table 4 in File SI & Figure 4E). Unigenes KA244417, 
KA275422 and KA275423 all matched AtMrB4, which is a key 
repressor of the phenylpropanoid pathway involved in the 
repression of lignin synthesis [42]. These three genes are 
upregulated in TW (Table 4 in File SI & Figure 4E); meanwhile, 
lignin content were both decreased in TW (Figure 1), which 
suggested that they may have similar functions as AtMYB4 in the 
repression of lignin synthesis. Unigene KA262304 matched 
AtMTB46, a key gene in the regulation of secondary wall 
biosynthesis in Arabidopsis [38]. As it is upregulated in TW, it 
may be involved in the formation of TW. The differences in AiYB 
expression between TW, OW and NW suggests that these MYB 
proteins may play positive or opposite roles in the process of TW 
formation. 

The Expression of Growth Hormone-related Genes 

Auxin is involved in acid-induced wall loosening and vascular 
dc\'clopmi'nt [43] and may be a key regulator of TW formation. A 
large number of auxin family genes were identified, including 
genes that encoded ARF (Auxin response factor), Aux/IAA 
protein, auxin efflux carrier protein, auxin influx transport protein, 
auxin response factor and auxin-induced protein (Table 4 in File 
SI). Most auxin family genes were differentially expressed among 
TW, NW and OW (Table 4 in File SI & Figure 4F), and were 
upregulated in TW compared in OW and NW, which implied that 
auxin-related genes are involved in the development of birch TW. 

Confirmation of the DGE Results by qRT-PCR Analysis 

To evaluate the validity of expression profiles detected by the 
transcriptome study, real-time RT-PCR analysis was performed. 
The real-time RT-PCR results were in agreement with the DGE 
results (Figure S3), which validated the reliability of the DGE 
results. 

Conclusions 

In this study, our results showed that the cellulose content was 
significantiy increased in TW compared with OW or NW, 
whereas the lignin content in TW was significantly decreased. 
Three transcriptome libraries were constructed from TW, OW 
and NW, and 9,684 significantly differentially expressed genes 
were identified. Among these, genes were identified that were 
involved in secondary cell wall structure and wood composition, 
including cellulose and lignin biosynthetic genes. The character- 
ization of gene expression in response to artificial bending on a 
genomic level will accelerate our understanding of the mechanisms 
that underlie xylem development in plants and will help us to 
identify genes involved in xylem development in birch. These 
findings may have potential applications for the improvement of 
wood properties in plants via genetic engineering. 

Supporting information 

Figure SI Result of GO analysis based on biological 
process (a), cellular components (b) and molecular 
function (c) of the three birch transcriptomes. 

(TIF) 



Figure S2 Pathway analysis of cellulose and lignin 
biosynthesis. Regulation of genes involved in cellulose (a) and 
lignin (b) synthesis during tension wood (TW) formation. The 
figure illustrates the sucrose, galactose, glucose and fructose 
metabolic pathways, which are related to cellulose synthesis and 
the shikimate, phenylpropanoid and monlignol biosynthetic 
pathways, which are related to lignin synthesis, according to 
http://www.genome.jp/kegg/pathway/map. Modifications of 
transcript abundance (Table 4 in File SI and Fig. 4) are indicated 
by the following colors: red, increase; blue, decrease; gray, not 
present in the three libraries; black, not affect(-d. Gc'n(;s that were 
significant for wood formation were identified according to signal 
strength and their relative abundance in the libraries. The 
following genes were included: cellulose synthase (CesA); sucrose 
synthase (SuSy); sucrose phosphate synthase (SPS); sucrose- 
phosphatase and sucrose phosphate phosphatase (SPP); UDP- 
glucose 4-epimerase (GALE); glucose- 1 -phosphate uridylyltrans- 
ferase; phosphoglucomutase; glucose-6-phosphate isomerase; 
hexokinase and fructokinase; beta-fructosidases; dehydroquinate 
dehydratase-shikimate dehydrogenase (DHQ;SDH), shikimate 
kinase (SK), ,'i-enolpyru\')'lshikimate-3-phosphate synthase 
(EPSPS); phenylalanine ammonialyase (PAL); 4-coumarate-coa 
ligase (4CL); cinnamyl alcohol dehydrogenase (COMT); caffeoyl- 
CoA O-methyltransferase (CCoAOMT); cinnamyl alcohol dehy- 
drogenase (CAD); Cinnamoyl-CoA reductase (CCR) and perox- 
idase (POD). 
(TIF) 

Figure S3 Confirmation of Solexa expression profiles by 
qRT-PCR analysis. NW was the control. AU ratios are log2 
transformed. 
(TIF) 

Table SI Differentially expressed genes identified by 
digital gene expression (DGE) analysis (>400 bp, FDR 
£0.0001, |log2Ratio| ^2). 

(XLS) 

File SI Contains Tables 1, 2, 3, and 4. Table 1. Primer 
sequences used for real time RT-PCR. Table 2. General 
characteristics of the three transcriptomes. Table 3. Results of 
GO analysis based on biological process, cellular components and 
molecular functions of the thrtx: birch transcriptomes. Table 4. 
Hierarchical cluster analysis of gene families that were differen- 
tially expressed in response to artificial bending. All gene families 
described in the Discussion are listed, including genes that were 
differentially expressed and those that were not. The differentially 
expressed genes were identified by a j log2Ratio j Si. 
(DOC) 
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